function StoreMarriage=solveMonteCarlo20191030(param,Glob)

alpha=param.alpha;
omega=param.omega;
sigma=param.sigma;
alpha1=alpha(:,:,1);
alpha2=alpha(:,:,2);
alpha3=alpha(:,:,3);

seed=Glob.seed;
dimS=Glob.dimS;
dimNation=Glob.dimNation;
dimProv=Glob.dimProv;
dimProv=length(omega);
dimNum=Glob.dimNum;

StoreMarriage=zeros(dimNation,dimNation,dimProv);

rng(seed,'twister');
typeM=randi(dimNation,dimS,dimNum,dimProv); 
typeF=randi(dimNation,dimS,dimNum,dimProv); 
typeMM=permute(typeM,[1 4 2 3]); 
typeMM=repmat(typeMM,[1 dimS 1 1]);
typeFF=permute(typeF,[4 1 2 3]); 
typeFF=repmat(typeFF,[dimS 1 1 1]); 
eps=randn(dimS,dimS,dimNum,dimProv)*sigma; 
love=exp(eps);
Alpha1=alpha1(typeMM(:,:,:,1)+(typeFF(:,:,:,1)-1)*dimNation); 
Alpha2=alpha2(typeMM(:,:,:,2)+(typeFF(:,:,:,2)-1)*dimNation); 
Alpha3=alpha3(typeMM(:,:,:,3)+(typeFF(:,:,:,3)-1)*dimNation); 
Alpha=cat(4,Alpha1,Alpha2,Alpha3);  
Omega=repmat(omega,[1 dimNum dimS dimS]);
Omega=permute(Omega,[3 4 2 1]);  
Surplus=love.*Alpha-2*Omega;  
parfor iter=1:(dimNum*dimProv)
    iprov=floor(iter/(dimNum+0.001))+1;
    inum=iter-dimNum*(iprov-1);

auxSurplus=Surplus(:,:,inum,iprov);    
matchresult=gurobSolve(max(auxSurplus,0),dimS);
match=reshape(matchresult.x,dimS,dimS);
[i1,OptWomen]=max(match,[],2);
ind=sub2ind([dimS dimS],(1:dimS)',OptWomen);
OptSurplus=auxSurplus(ind); 
SinglesM=OptSurplus<0; 
Pred_Marriage=zeros(dimNation,dimNation);  
M=[typeM(:,iter) typeF(OptWomen,iter)];
q=sub2ind([dimNation dimNation],M(not(SinglesM),1),M(not(SinglesM),2));
tab=accumarray(q,1);
Pred_Marriage(1:length(tab))=tab;
aux=zeros(dimNation,dimNation,dimProv); 
aux(:,:,iprov)=Pred_Marriage;
StoreMarriage=StoreMarriage+aux;
end
StoreMarriage=StoreMarriage/dimS/dimNum;
